.libPaths("/dora/nobackup/yuec/R")
library(SingleCellExperiment)
library(ggplot2)
library(scFeatures)
library(ClassifyR)
library(lisaClust)
library(ggthemes)
library(spicyR)
library(dplyr)
library(limma)
library(plotly)
library(scattermore)
library(tidyr)
library(survival)
library(survminer)
library(spatstat.core)
library(spatstat.geom)
library(scater)
library(scran)
library(SPOTlight)
library(reshape)
theme_set(theme_classic())

1 Overview

As single cell technology advances, the recent development of spatial omics allows us to examine the spatial organisation of cells within tissue in their native environment. This workshop will discuss the challenges and analytical focus associated with disease outcome prediction using multi-condition multi-sample spatial dataset. We will also talk about general analytic strategies and the critical thinking questions that arise in the workflow.



1.0.1 Preparation and assumed knowledge

1.0.2 Learning objectives

  • Understand and visualise spatial omics dataset
  • Explore various strategies for disease outcome prediction using spatial omics data
  • Understand the transformation from cell level features to patient level features
  • Generate patient representations from gene expression matrix
  • Understand the characteristics of good classification models
  • Perform disease outcome prediction using the feature representation and robust classification framework



2 Background of the dataset

In this demo, we look at a Visium dataset taken from Kuppe, C., Ramirez Flores, R. O., Li, Z., Hayat, S., Levinson, R. T., Liao, X., … & Kramann, R. (2022). Spatial multi-omic map of human myocardial infarction. Nature, 608(7924), 766-777.

Visium captures spatial information, creating images that display the distribution of different cell types and their associated gene expression patterns in the tissue.

In this dataset, the authors quantified the expression of >10000 genes in control and in patients with myocardial infarction. In this demo, we examine patients defined to be myogenic group and ischaemic group. Myogenic group is defined by sample taken from control and unaffected remote zone, ischaemic group is defined by sample taken from ischaemic zone.

3 Initial exploration and visualisation

3.1 Exploration 1: How complex is my data ?

Examine the data objects:
- The dataset contains 11,681 genes and 19,000 cells (after subsampling).
- The outcome is 11 myogenic samples (4 from control, 7 from remote zone) and 8 ischaemic samples.

data_sce <- readRDS("~/data/small_data.rds")
 
data_sce
## class: SpatialExperiment 
## dim: 11681 19000 
## metadata(3240): X_approximate_distribution X_normalization ...
##   schema_version title
## assays(2): X logcounts
## rownames(11681): LINC01409 LINC01128 ... VAMP7 AC007325.4
## rowData names(1): features
## colnames(19000): ACAGAACTGAGAACAA-1 TTGCCCTGATCACGGG-1 ...
##   CTGACATAGAAATAGA-1.16 ATTAGGCGATGCTTTC-1.24
## colData names(33): n_counts n_genes ... number_cells condition
## reducedDimNames(3): X_pca X_spatial X_umap
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(0) :
## imgData names(1): sample_id
print("expression matrix is stored in genes by cells matrix")
## [1] "expression matrix is stored in genes by cells matrix"
logcounts(data_sce)[1:7, 1:7]
## 7 x 7 sparse Matrix of class "dgCMatrix"
##           ACAGAACTGAGAACAA-1 TTGCCCTGATCACGGG-1 CATATAGGTACAGTCA-1
## LINC01409           .                  .                  .       
## LINC01128           1.431438           .                  1.351835
## NOC2L               .                  .                  .       
## KLHL17              .                  .                  .       
## HES4                .                  .                  .       
## ISG15               .                  1.279462           .       
## AGRN                .                  .                  .       
##           TGCGAGAATATTACCC-1 TATTCAATTCTAATCC-1 GTCCTATTGTTGTGGT-1
## LINC01409           .                  .                  .       
## LINC01128           .                  .                  .       
## NOC2L               .                  .                  .       
## KLHL17              .                  .                  .       
## HES4                1.146986           1.389528           .       
## ISG15               .                  .                  1.454652
## AGRN                .                  .                  .       
##           TAACCTACCGTCCGAG-1
## LINC01409           .       
## LINC01128           .       
## NOC2L               .       
## KLHL17              .       
## HES4                .       
## ISG15               1.475897
## AGRN                .
print("the object stores meta data (such as patient outcome information) about each cell")
## [1] "the object stores meta data (such as patient outcome information) about each cell"
DT::datatable( data.frame(colData(data_sce))[1:5, ], options = list(scrollX = TRUE))

3.2 Exploration 2: How to visualise my data ?

Typically in a spatial data, we perform dimension reduction to reduce and project the high dimensional cell by gene matrix on to 2D space. This allows us to visualise various things of interest, such as distribution of cell types and disease outcomes.

Visium is a spot-based technology, meaning that each spot capture 1-10 cells. Therefore each spot may represent a mixture of cells from different cell types. In this dataset, the author used cell2location which is a spatial deconvolution tool to predict the cell type probability or composition of each of the spot. There are also many other deconvolution tools available, eg, CARD published in Ma, Y., & Zhou, X. (2022). Spatially informed cell-type deconvolution for spatial transcriptomics. Nature biotechnology, 40(9), 1349-1359.

For demonstration purpose, in the plot below we visualise each spot using the cell type with maximum probability.

data_sce <- runUMAP(data_sce, exprs_values = "logcounts", scale=TRUE, min_dist = 0.3)
a <- plotUMAP(data_sce, colour_by = "celltype")
b <- plotUMAP(data_sce, colour_by = "condition")
c <- plotUMAP(data_sce, colour_by = "sample")
ggarrange(plotlist = list(a,b,c))

Interactive Q&A:

  • Q1: Is there patient batch effect?
  • Q2: Are the myogenic group and ischaemic group easy or difficult to distinguish?
  • Any interesting patterns from the plot?

3.3 Exploration 3: Is there a spatial structure in my data?

The advantage with spatial omics is that we can examine the organisation of the cell types as it occurs on the tissue slide. Here we visualise one of the slides from a selected patient. As an optional exercise, you can:

  • permute the cell type label.
  • permute the spatial coordinate.

to give a sense of what is random ordering.

celltype <- c("Adipocyte" , "Cardiomyocyte" , 
              "Endothelial" ,  "Fibroblast" , "Lymphoid" ,    
              "Mast" ,   "Myeloid" , 
              "Neuronal" , "Pericyte" , "Cycling.cells", "vSMCs" )

tableau_palette <- scale_colour_tableau() 
color_codes <-  c( tableau_palette$palette(10) , "lightgrey")
names(color_codes) <- celltype


 
one_sample_data <- data_sce[, data_sce$sample  ==  "IZ_P9_cond_Ischaemic"]
one_sample  <- colData(one_sample_data)
one_sample <- data.frame(one_sample)
a <- ggplot(one_sample, aes(x = spatial_x , y = spatial_y, colour = celltype)) + geom_point(alpha=0.7) +  scale_colour_manual(values = color_codes) + ggtitle("Original slide")
print("Optional: Permute the cell type labels")
## [1] "Optional: Permute the cell type labels"
one_sample$celltype_permute <- sample(one_sample$celltype)
b <- ggplot(one_sample, aes(x = spatial_x , y =  spatial_y, colour =celltype_permute)) + geom_point(alpha=0.7) +  scale_colour_manual(values = color_codes)  + ggtitle("Permute the cell type label")
print("Optional: Permute the spatial coordinate")
## [1] "Optional: Permute the spatial coordinate"
one_sample$spatial_x_permute <- sample(one_sample$spatial_x)
one_sample$spatial_y_permute <- sample(one_sample$spatial_y)

c <- ggplot(one_sample, aes(x = spatial_x_permute , y = spatial_y_permute, colour = celltype)) + geom_point(alpha=0.7) +  scale_colour_manual(values = color_codes) + ggtitle("Permute the X, Y coordinate")


ggarrange(plotlist = list(a,b,c), ncol =3 )

Instead of plotting the cell type with maximum probability, we can also visualise the cell type composition of each spot using pie chart.

x <- data.frame( imagecol = one_sample$spatial_y, 
                 imagerow = one_sample$spatial_x)
 
rownames(x) <- paste0("Spot", 1:nrow(x))

y <- data.frame( colData(one_sample_data)[ ,celltype  ] )
 
rownames(y) <- paste0("Spot", 1:nrow(y))

plotSpatialScatterpie(x = x, y = y , pie_scale = 0.7)+theme_classic() +  scale_fill_manual(values = color_codes) + ylab("spatial_y") + xlab("spatial_y")

Interactive Q&A:

Q3: Is there a structure in the data or is the cell type randomly distribution?

4 Describing tissue microenvrionments and cellular neighbourhoods

4.1 Do cell type co-localise in specfic regions?

We begin by examine how can we identify and visualise regions of tissue where spatial associations between cell-types is similar? There are many packages that perform this task andhere we use the lisaClust function [https://www.bioconductor.org/packages/devel/bioc/html/lisaClust.html] that is based on “local L-function” to spatially clusters cells into different regions with similar cell type composition.

set.seed(51773)
 

BPPARAM <- simpleSeg:::generateBPParam(2)

# Cluster cells into spatial regions with similar composition.
data_sce  <- lisaClust(
  data_sce ,
  k = 5,
  Rs = c(20, 50, 100),
  sigma = 50,
  spatialCoords = c("spatial_x", "spatial_y"),
  cellType = "celltype",
  imageID = "sample" ,
  regionName = "region",
  BPPARAM = BPPARAM
)

4.2 Which regions appear to be different between the ischaemic and myogenic patients?

4.2.1 Visualise regions on individual level

Using the slide we previously shown, visualise the region output.

metadata <- colData(data_sce)
metadata <- metadata[ metadata$sample == "IZ_P9_cond_Ischaemic",  ]
metadata <- data.frame(metadata)
metadata$celltype  <- as.character(metadata$celltype )

plotlist <- list()
plotlist_celltype <- list()
thisregion  <-  unique(metadata$region)[1]


tableau_palette <- scale_colour_tableau() 
color_codes <- tableau_palette$palette( 10 )
color_codes <- c(color_codes, "darkgrey" , "grey90") 

names(color_codes) <- c( celltype ,  "other regions")

for ( thisregion in sort(unique(metadata$region))){
        
        selected_region_index <-  metadata$region == thisregion
        
        metadata$selected_region <-  "other regions"
        metadata$selected_region[selected_region_index] <- "selected region"
        
        metadata$selected_celltype <- metadata$celltype
        metadata$selected_celltype[!selected_region_index ] <-   "other regions"
        
        # metadata$celltype <- factor(metadata$celltype, levels = c(unique(metadata$celltype), "other regions"))

       p <- ggplot(metadata, aes(x = spatial_x , y = spatial_y, colour = selected_region  )) + 
         geom_point( alpha = 0.8 ) + ggtitle(thisregion) + scale_colour_manual(values = c("grey" , "red"))
         
       
    
       p2 <-  ggplot(metadata, aes(x = spatial_x , y = spatial_y, colour =  selected_celltype )) +
         geom_point(alpha = 0.8 ) + ggtitle(thisregion) + scale_colour_manual(values =  color_codes)
       
      plotlist [[thisregion ]] <- p
       
      plotlist_celltype [[thisregion ]] <- p2
}

ggarrange(plotlist = plotlist , ncol = 5, nrow = 1 , common.legend = T )

ggarrange(plotlist = plotlist_celltype , ncol = 5, nrow = 1 , common.legend = T )

4.2.2 Visualise regions across patients

We can better interpret the region output by summarising the proportion of each cell type in a region across the patients.

Looking at the composition of cell types in each region, comparing between ischaemic and myogenic.

df <- data.frame(colData( data_sce))
 

df_plot <- NULL
for ( thispatient in unique(df$sample)){
  this_df <- df[df$sample == thispatient, ]
  temp_df <-   table( this_df$region, this_df$celltype )
  temp_df <-  temp_df / rowSums(temp_df)
   temp_df <- data.frame(  temp_df)
  temp_df$patient <-  thispatient
  temp_df$reponse <- unique( this_df$condition)
  df_plot <- rbind(df_plot, temp_df)
}

df_plot <- df_plot %>% dplyr::group_by( Var1 , Var2, reponse) %>% 
  summarise(mean_proportion = mean(Freq))
  
 
ggplot(df_plot, aes(y = Var2, x = Var1 ,colour =mean_proportion  , size = mean_proportion ))+  geom_point() + 
  facet_grid(~reponse, scales = "free", space = "free" ) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + 
  xlab("Region" ) + ylab("Celltype") + scale_colour_viridis_c()

Interactive Q&A:

Q4: Which regions would you focus on next ?

df <- data.frame(colData( data_sce))

df <- df %>% dplyr::group_by(sample , condition , region) %>%
  count(celltype) %>%
  mutate(proportion = n / sum(n))


ggplot(df, aes(y = proportion, x = sample , fill = celltype))+ geom_col()+facet_grid(~region+condition, scales = "free", space = "free" ) + scale_fill_manual(values = c(color_codes))  +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

Interactive Q&A:

Q5: Does your conclusion change after looking at a different plot

4.3 Further exploration by visualise selected regions

Let’s focus on region 2 and visualise the boxplot of cell type distribution.

df <- data.frame(colData( data_sce))



df_plot <- NULL
for ( thispatient in unique(df$sample )){
  this_df <- df[df$sample == thispatient, ]
  temp_df <-   table( this_df$region, this_df$celltype)
  temp_df <-  temp_df / rowSums(temp_df)
  temp_df <- data.frame(  temp_df)
  temp_df$patient <-  thispatient
  temp_df$reponse <- unique( this_df$condition )
  df_plot <- rbind(df_plot, temp_df)
}

df_plot_region_1 <- df_plot[df_plot$Var1 == "region_2", ]
 
a <- ggplot(df_plot_region_1, aes(x =  Var2,  y = Freq, colour = reponse)) +
  geom_boxplot()+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + 
  ylab("Proportion") + xlab("Cell type")+ ggtitle("Region 2") + ylim(0,1)

a

Discussion:

Comparing the cell type composition in the region, what can you tell about the regions?

5 How do we generate molecular representation for each individual ?

In this demo, we use scFeatures to generate molecular representation for each patient. The molecular representation is interpretable and hence facilitates downstream analysis of the patient. Overall, scFeatures generates features across six categories representing different molecular views of cellular characteristics. These include: - i) cell type proportions - ii) cell type specific gene expressions - iii) cell type specific pathway expressions - iv) cell type specific cell-cell interaction (CCI) scores - v) overall aggregated gene expressions - vi) spatial metrics The different types of features constructed enable a more comprehensive multi-view understanding of each patient from a matrix of proteins x cells.

print("number of cells in each sample")
## [1] "number of cells in each sample"
metadata <- data.frame( table(data_sce$sample) )
colnames(metadata)[1] <- "Patient"
DT::datatable(metadata , options = list(pageLength = 5), width = "400px")
print("number of cells in each celltype")
## [1] "number of cells in each celltype"
metadata <- data.frame( table(data_sce$celltype) )
colnames(metadata)[1] <- "Region specific cell type"
DT::datatable(metadata , options = list(pageLength = 5), width = "400px")

Discussion:

Is there any sample or cell types you would like to remove from the data?

5.1 How to create one molecular representations for an individual ?

All the feature types can be generated in one line of code. This runs the function using default settings for all parameters, for more information, type ?scFeatures.

#scFeatures requires the following columns 
#  data, sample, x_cord and y_cord
data <- logcounts( data_sce )
sample <- data_sce$sample
spatialCoords <- list(data_sce$spatial_x, data_sce$spatial_y)


# for spatial transcriptomics, we also need to specify the cell type prediction of each spot 
prediction.scores <- colData(data_sce)[, celltype]
prediction.scores <- as.matrix( t(prediction.scores)  )

 
scfeatures_result <- scFeatures(data = data ,  sample = sample, 
                spatialCoords = spatialCoords , spotProbability = prediction.scores, 
                type ="spatial_t", ncores = 20)

5.2 Visualising and exploring scFeatures output

We have generated a total of 13 feature types and stored them in a list. All generated feature types are stored in a matrix of samples by features. For example, the first list element contains the feature type “proportion_raw”, which contains the cell type proportion features for each patient sample. We could print out the first 5 columns and first 5 rows of the first element to see.

scfeatures_result <- readRDS("~/data/scfeatures_result_matrix.rds")
 
# we have generated a total of 13 feature types
names(scfeatures_result )
##  [1] "proportion_raw"       "proportion_logit"     "proportion_ratio"    
##  [4] "gene_mean_celltype"   "pathway_mean"         "gene_mean_bulk"      
##  [7] "gene_prop_bulk"       "gene_cor_bulk"        "L_stats"             
## [10] "celltype_interaction" "morans_I"             "nn_correlation"
# each row is a sample, each column is a feature 

meta_table <- data.frame(round ( scfeatures_result [[1]][1:5, 1:5] , 2))
DT::datatable(meta_table , options = list(scrollX = TRUE))

Once the features are generated, you may wish to visually explore the features.

Here we plot a volcano plot and a dotplot for the cell type specific expression feature.

gene_mean_bulk <- scfeatures_result$gene_mean_bulk
# this transposes the data
# in bioinformatics conversion, features are stored in rows 
# in statistics convention, features are stored in columns
gene_mean_bulk <- t(gene_mean_bulk)
      
 
condition  <- unlist( lapply( strsplit( colnames(gene_mean_bulk), "_cond_"), `[`, 2))
condition <- data.frame(condition = condition )
design <- model.matrix(~condition, data = condition)
fit <- lmFit(gene_mean_bulk, design)
fit <- eBayes(fit)
tT <- topTable(fit, n = Inf)
tT$gene <- rownames(tT)
p <- ggplot( tT , aes(logFC,-log10(P.Value) , text = gene ) )+
      geom_point(aes(colour=-log10(P.Value)), alpha=1/3, size=1) +
      scale_colour_gradient(low="blue",high="red")+
      xlab("log2 fold change") + ylab("-log10 p-value")  
 
p

tT <- tT[ order(tT$logFC, decreasing = T), ]
tT <- tT[1:20, ]
ggplot( tT , aes( y = reorder(gene, logFC) , x = logFC  ) )+
      geom_point(aes(colour=-log10(P.Value)), alpha=1/3, size=4) +
      scale_colour_gradient(low="blue",high="red")+
      xlab("logFC") + ylab("region specific cel type specfic features" )

Interactive Q&A:

Q6: Which figure do you prefer? The volcano plot or the dotplot?

We can also visualise the distribution of the features across the patients. This helps to examine whether there exists any batch effect due to patient effect.

plot_boxplot <- function( feature ){
  
  data <- t(feature)

  patient <- unlist( lapply( strsplit( colnames(data), "_cond_"), `[`, 1) ) 
  condition  <- unlist( lapply( strsplit( colnames(data), "_cond_"), `[`, 2))
  condition <- data.frame(condition = condition )
      
  design <- model.matrix(~condition, data = condition)
  fit <- lmFit(data, design)
  fit <- eBayes(fit)
  tT <- topTable(fit, n = Inf) 
  
  
  # selecting the top 10 DE features 
  top_gene <- rownames( tT)[1:10 ]
  rownames( condition) <- colnames(data)
   
  data_plot <-  data[top_gene, ]   
  
  data_plot <- melt(data_plot )
  data_plot$condition <- unlist( lapply( strsplit(  as.character( data_plot$X2), "_cond_"), `[`, 2))
  
  
  p <- ggplot(data_plot, aes( x = X1, y = value , colour = condition)) + 
    geom_boxplot() + geom_point() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
  
   
  return(p)
  
}
  


plot_barplot <- function(data , dodge=F  ){
  

  data$patient <- unlist( lapply( strsplit( rownames(data ), "_cond_"), `[`, 1))
  data$condition <- unlist( lapply( strsplit( rownames(data ), "_cond_"), `[`, 2))
  
  data <- as.data.frame( melt(data, id=c("patient", "condition")) )
 
 
  p <-   ggplot(data , aes( x = patient , y = value , fill = variable) ) +   
    geom_bar(stat="identity"   ) + facet_wrap(~condition, scale="free") + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))  + 
    scale_fill_manual(values = color_codes)
 
  
 
 return (p)
  
 
 
}
feature <- scfeatures_result$proportion_raw
 
a <- plot_barplot(feature ) + ggtitle("Proportion raw feature")


feature <- scfeatures_result$gene_mean_celltype
feature   <- feature[  , grep("Cardio", colnames(feature)) ]
b <- plot_boxplot(feature) + ggtitle("Gene mean celltype feature")

feature <- scfeatures_result$pathway_mean
c <- plot_boxplot(feature) + ggtitle("Pathway mean feature")

feature <- scfeatures_result$gene_mean_bulk
d <- plot_boxplot(feature) + ggtitle("Gene mean bulk feature")

feature <- scfeatures_result$nn_correlation
e <- plot_boxplot(feature) + ggtitle("Nearest neighbour correlation feature")


ggarrange(plotlist = list(a, b,c,d,e) , ncol = 2, nrow = 3 )

To accommodate for easier interpretation of the features, scFeatures contains a function run_association_study_report that enables the user to readily visualise and explore all generated features with one line of code.

# specify a folder to store the html report. Here we store it in the current working directory. 
output_folder <-  getwd()
run_association_study_report(scfeatures_result, output_folder )

5.3 Are the generated features sensible ?

Interactive Q&A:

Using the HTML, we can look at some of the critical thinking questions that a researcher would ask about the generated features. These questions are exploratory and there is no right or wrong answer.

Q7: Do the generated features look reasonable?
Which cell type(s) would you like to focus on at your next stage of analysis?
Which feature type(s) would you like to focus on at your next stage of analysis?
Q8: Are the conditions in your data relatively easy or difficult to distinguish?

6 Can we classify or discrimiante between the ischaemic and myogenic patients?

6.1 Building classification model

Recall in the previous section that we have stored the 13 feature types matrix in a list. Instead of manually retrieving each matrix from the list to build separate models, classifyR can directly take a list of matrices as an input and run repeated cross-validation model on each matrix individually.

Below, we run 5 repeats of 3 folds cross-validation.

outcome = scfeatures_result[[1]] %>% rownames %>% strsplit("_cond_") %>% sapply(function(x) x[2])
table(outcome)

### generate classfyr result 
classifyr_result <- crossValidate(scfeatures_result,
                                 outcome, 
                                 classifier = "kNN",
                                 nFolds = 3, 
                                 nRepeats = 5, 
                                 nCores = 5  )

6.2 Visualising the classification performance

To examine the classification model performance, we first need to specify a metric to calculate. Here, we calculate the balanced accuracy.

classifyr_result <-  readRDS("~/data/classifyr_result.rds")

classifyr_result <- lapply(classifyr_result, 
                           function(x) calcCVperformance(x, performanceType = "Balanced Accuracy"))

Format the output and visualise the accuracy using boxplots.

level_order <- names(scfeatures_result)

p  <- performancePlot(classifyr_result) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + 
  scale_x_discrete(limits = level_order)  

p

7 Appendix

7.1 Explanation of spatial features

  • L function:

The L function is a spatial statistic used to assess the spatial distribution of cell types. It assess the significance of cell-cell interactions, by calculating the density of a cell type with other cell types within a certain radius. High values indicate spatial association, low values indicate spatial avoidance.

tableau_palette <- scale_colour_tableau() 
color_codes <- tableau_palette$palette( 10 )
color_codes <- c(color_codes, "darkgrey" , "grey90") 

names(color_codes) <- c( celltype ,  "other regions")
 


one_sample_data  <- data_sce[ , data_sce$sample == "IZ_P10_cond_Ischaemic"  ]

one_sample <- data.frame( colData(one_sample_data) )
index <-  one_sample$celltype  %in% c("Cardiomyocyte", "Myeloid")
one_sample$celltype <- as.character(one_sample$celltype)
one_sample$celltype[!index] <- "others"
a <- ggplot( one_sample, aes(x = spatial_x , y = spatial_y, colour = celltype )) + geom_point()  + scale_colour_manual(values = c("steelblue", "coral", "grey"))  + ggtitle( "Patient P10 - high L value with \n cardiomyocyte interacting myeloid")
 
one_sample <- data.frame( colData(one_sample_data) )
index <-  one_sample$celltype  %in% c("Cardiomyocyte", "Cycling.cells")
one_sample$celltype <- as.character(one_sample$celltype)
one_sample$celltype[!index] <- "others"
b <- ggplot( one_sample, aes(x = spatial_x , y = spatial_y, colour = celltype )) + geom_point() + scale_colour_manual(values = c("steelblue", "coral", "grey"))  + ggtitle( "Patient 10 - low L value with  \n cardiomyocyte interacting cycling cells")
 
ggarrange(plotlist = list( a , b))

  • Cell type interaction composition:

We calculate the nearest neighbours of each cell and then calculate the pairs of cell type based on the nearest neighbour. This allow us to summarised it into a cell type interaction composition.

tableau_palette <- scale_colour_tableau() 
color_codes <- tableau_palette$palette( 10 )
color_codes <- c(color_codes, "darkgrey" , "grey90") 

names(color_codes) <- c( unique(data_sce$celltype)  )
 

tableau_palette <- scale_colour_tableau() 
color_codes <-  c( tableau_palette$palette(10) , "lightgrey")
names(color_codes) <- celltype



one_sample  <- data_sce[ , data_sce$sample == "IZ_P10_cond_Ischaemic"  ]
one_sample <- data.frame( colData(one_sample) )
 
a <- ggplot( one_sample, aes(x = spatial_x , y = spatial_y, colour = celltype )) + geom_point( alpha = 0.8)  + scale_colour_manual(values = color_codes)  + ggtitle( "Patient P10")


feature  <- scfeatures_result$celltype_interaction
to_plot <- data.frame( t( feature["IZ_P10_cond_Ischaemic" , ])  )
to_plot$feature <- rownames(to_plot) 
colnames(to_plot)[2] <- "celltype interaction composition"
 to_plot <- to_plot[ order(to_plot$IZ_P10_cond_Ischaemic , decreasing = T), ]
 to_plot <-  to_plot[1:5 , ]
 # to_plot <- to_plot[ to_plot$IZ_P10_cond_Ischaemic  > 0.03 , ] 
b <- ggplot(to_plot, aes( x =  reorder(`celltype interaction composition`, IZ_P10_cond_Ischaemic) ,  y = IZ_P10_cond_Ischaemic  ,fill = `celltype interaction composition`)) + geom_bar(stat="identity", fill = "steelblue" ) + ylab("Cell type interactions value")  + xlab("Major cell type interactions") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) 

ggarrange(plotlist = list( a , b))

  • Moran’s I:

Moran’s I is a spatial autocorrelation statistic based on both location and values. It quantifies whether similar values tend to occur near each other or dispersed.

high  <- data_sce["PDGFD", data_sce$sample == "IZ_P16_cond_Ischaemic"  ]
high_meta <- data.frame( colData(high) ) 
high_meta$expression <- as.vector(logcounts( high)) 

low  <- data_sce["KIF22"  , data_sce$sample == "IZ_P16_cond_Ischaemic" ]
low_meta <- data.frame( colData(low) )
low_meta$expression <- as.vector(logcounts(low))


a <- ggplot(low_meta, aes(x = spatial_x , y = spatial_y, colour =expression)) + geom_point(alpha=0.5) + scale_colour_viridis_c() + ggtitle("Patient P16 - low Moran's I in KIF22")

b <- ggplot(high_meta, aes(x = spatial_x , y = spatial_y, colour =expression)) + geom_point(alpha=0.5) + scale_colour_viridis_c() + ggtitle("Patient P16 - high Moran's I in PDGFD")

ggarrange(plotlist = list( a , b))

  • Nearest Neighbor Correlation:

This metric measures the correlation of proteins/genes between a cell and its nearest neighbor cell.

We see from both prediction and survival analysis that the feature type “nn correlation” (nearest neighbour correlation) perform the best.

Here we pick the gene “COL1A1” as an example to illustrate the concept.

 plot_nncorrelation <- function(thissample , thisprotein){
   
       sample_name <- thissample
       thissample <- data_sce[, data_sce$sample ==     sample_name]
    
      
      exprsMat <- logcounts(thissample)
     
    
    cell_points_cts <- spatstat.geom::ppp(
            x = as.numeric(thissample$spatial_x), y = as.numeric(thissample$spatial_y),
            check = FALSE,
            xrange = c(
                min(as.numeric(thissample$spatial_x)),
                max(as.numeric(thissample$spatial_x))
            ),
            yrange = c(
                min(as.numeric(thissample$spatial_y)),
                max(as.numeric(thissample$spatial_y))
            ),
            marks = t(as.matrix(exprsMat))
        )
    
     value <-  spatstat.explore::nncorr(cell_points_cts)["correlation", ]
      value <-  value[  thisprotein]
     
    # Find the indices of the two nearest neighbors for each cell
    nn_indices <- nnwhich(cell_points_cts, k = 1)
    
    protein <-  thisprotein
    df <- data.frame(thiscell_exprs  = exprsMat[protein, ] , exprs =  exprsMat[protein,nn_indices ])
    
   p <-  ggplot(df, aes( x =thiscell_exprs ,  y = exprs , colour =  exprs  )) +
      geom_point(alpha = 0.3) + ggtitle(paste0( "Patient ", sample_name ,  " nn_corr = " ,  round(value, 2)  )) + scale_colour_viridis_c() + xlab("This cell expression") + ylab("Neighbouring cell expression")
   
   return (p ) 

}

    
p1 <- plot_nncorrelation( "IZ_P9_cond_Ischaemic" ,  "COL1A1" )
p2 <- plot_nncorrelation( "control_P1_cond_Myogenic"  ,  "COL1A1" )
ggarrange(plotlist = list( p1 , p2))

7.2 Session info

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 12 (bookworm)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.21.so;  LAPACK version 3.11.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: Australia/Sydney
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] reshape_0.8.9               SPOTlight_1.5.2            
##  [3] scran_1.28.1                scater_1.28.0              
##  [5] scuttle_1.10.1              spatstat.core_2.3-0        
##  [7] rpart_4.1.19                nlme_3.1-162               
##  [9] spatstat.geom_3.2-1         spatstat.data_3.0-1        
## [11] survminer_0.4.9             ggpubr_0.6.0               
## [13] tidyr_1.3.0                 scattermore_0.8            
## [15] plotly_4.10.2               limma_3.56.2               
## [17] dplyr_1.1.2                 spicyR_1.12.0              
## [19] ggthemes_4.2.4              lisaClust_1.9.2            
## [21] ClassifyR_3.4.8             survival_3.5-5             
## [23] BiocParallel_1.34.2         MultiAssayExperiment_1.26.0
## [25] generics_0.1.3              scFeatures_0.99.27         
## [27] ggplot2_3.4.2               SingleCellExperiment_1.22.0
## [29] SummarizedExperiment_1.30.2 Biobase_2.60.0             
## [31] GenomicRanges_1.52.0        GenomeInfoDb_1.36.1        
## [33] IRanges_2.34.1              S4Vectors_0.38.1           
## [35] BiocGenerics_0.46.0         MatrixGenerics_1.12.2      
## [37] matrixStats_1.0.0          
## 
## loaded via a namespace (and not attached):
##   [1] SpatialExperiment_1.11.2   R.methodsS3_1.8.2         
##   [3] GSEABase_1.62.0            progress_1.2.2            
##   [5] tiff_0.1-11                EnsDb.Mmusculus.v79_2.99.0
##   [7] DT_0.28                    goftest_1.2-3             
##   [9] Biostrings_2.68.1          HDF5Array_1.28.1          
##  [11] vctrs_0.6.1                spatstat.random_3.1-5     
##  [13] digest_0.6.33              png_0.1-8                 
##  [15] shape_1.4.6                EnsDb.Hsapiens.v79_2.99.0 
##  [17] registry_0.5-1             ggrepel_0.9.3             
##  [19] deldir_1.0-9               parallelly_1.36.0         
##  [21] magick_2.7.4               MASS_7.3-59               
##  [23] reshape2_1.4.4             httpuv_1.6.11             
##  [25] foreach_1.5.2              withr_2.5.0               
##  [27] ggfun_0.1.1                xfun_0.39                 
##  [29] ellipsis_0.3.2             cytomapper_1.12.0         
##  [31] memoise_2.0.1              proxyC_0.3.3              
##  [33] ggbeeswarm_0.7.2           systemfonts_1.0.4         
##  [35] zoo_1.8-12                 GlobalOptions_0.1.2       
##  [37] gtools_3.9.4               SingleCellSignalR_1.12.0  
##  [39] pbapply_1.7-2              R.oo_1.25.0               
##  [41] prettyunits_1.1.1          KEGGREST_1.40.0           
##  [43] promises_1.2.0.1           httr_1.4.6                
##  [45] rstatix_0.7.2              restfulr_0.0.15           
##  [47] globals_0.16.2             fitdistrplus_1.1-11       
##  [49] rhdf5filters_1.12.1        rhdf5_2.44.0              
##  [51] rstudioapi_0.14            miniUI_0.1.1.1            
##  [53] concaveman_1.1.0           babelgene_22.9            
##  [55] curl_5.0.2                 zlibbioc_1.46.0           
##  [57] ScaledMatrix_1.8.1         polyclip_1.10-4           
##  [59] GenomeInfoDbData_1.2.10    fftwtools_0.9-11          
##  [61] doParallel_1.0.17          xtable_1.8-4              
##  [63] stringr_1.5.0              evaluate_0.21             
##  [65] S4Arrays_1.0.6             BiocFileCache_2.8.0       
##  [67] hms_1.1.3                  irlba_2.3.5.1             
##  [69] colorspace_2.1-0           filelock_1.0.2            
##  [71] ROCR_1.0-11                reticulate_1.30           
##  [73] magrittr_2.0.3             lmtest_0.9-40             
##  [75] viridis_0.6.3              later_1.3.1               
##  [77] lattice_0.21-8             NMF_0.26                  
##  [79] future.apply_1.11.0        XML_3.99-0.14             
##  [81] cowplot_1.1.1              RcppAnnoy_0.0.21          
##  [83] svgPanZoom_0.3.4           class_7.3-22              
##  [85] pillar_1.9.0               simpleSeg_1.2.0           
##  [87] EBImage_4.42.0             iterators_1.0.14          
##  [89] gridBase_0.4-7             caTools_1.18.2            
##  [91] compiler_4.3.1             beachmat_2.16.0           
##  [93] stringi_1.7.12             tensor_1.5                
##  [95] minqa_1.2.5                GenomicAlignments_1.36.0  
##  [97] plyr_1.8.8                 msigdbr_7.5.1             
##  [99] crayon_1.5.2               abind_1.4-5               
## [101] BiocIO_1.10.0              locfit_1.5-9.8            
## [103] sp_2.0-0                   terra_1.7-39              
## [105] bit_4.0.5                  codetools_0.2-19          
## [107] BiocSingular_1.16.0        crosstalk_1.2.0           
## [109] bslib_0.5.0                multtest_2.56.0           
## [111] mime_0.12                  splines_4.3.1             
## [113] circlize_0.4.15            Rcpp_1.0.10               
## [115] dbplyr_2.3.2               sparseMatrixStats_1.12.2  
## [117] knitr_1.43                 blob_1.2.4                
## [119] utf8_1.2.3                 AnnotationFilter_1.24.0   
## [121] lme4_1.1-34                nnls_1.4                  
## [123] listenv_0.9.0              DelayedMatrixStats_1.22.1 
## [125] GSVA_1.48.2                ggsignif_0.6.4            
## [127] tibble_3.2.1               Matrix_1.5-4.1            
## [129] scam_1.2-14                statmod_1.5.0             
## [131] svglite_2.1.1              tweenr_2.0.2              
## [133] pkgconfig_2.0.3            pheatmap_1.0.12           
## [135] tools_4.3.1                cachem_1.0.8              
## [137] RSQLite_2.3.1              viridisLite_0.4.2         
## [139] DBI_1.1.3                  numDeriv_2016.8-1.1       
## [141] fastmap_1.1.1              rmarkdown_2.23            
## [143] scales_1.2.1               grid_4.3.1                
## [145] ica_1.0-3                  Seurat_4.3.0.1            
## [147] shinydashboard_0.7.2       Rsamtools_2.16.0          
## [149] broom_1.0.5                sass_0.4.6                
## [151] patchwork_1.1.2            BiocManager_1.30.22       
## [153] graph_1.78.0               carData_3.0-5             
## [155] RANN_2.6.1                 farver_2.1.1              
## [157] scatterpie_0.2.1           mgcv_1.8-42               
## [159] yaml_2.3.7                 rtracklayer_1.60.0        
## [161] cli_3.6.1                  purrr_1.0.1               
## [163] leiden_0.4.3               lifecycle_1.0.3           
## [165] uwot_0.1.16                bluster_1.10.0            
## [167] backports_1.4.1            annotate_1.78.0           
## [169] gtable_0.3.4               rjson_0.2.21              
## [171] ggridges_0.5.4             progressr_0.13.0          
## [173] parallel_4.3.1             ape_5.7-1                 
## [175] jsonlite_1.8.7             edgeR_3.42.4              
## [177] bitops_1.0-7               bit64_4.0.5               
## [179] Rtsne_0.16                 spatstat.utils_3.0-3      
## [181] BiocNeighbors_1.18.0       SeuratObject_4.1.3        
## [183] RcppParallel_5.1.7         highr_0.10                
## [185] jquerylib_0.1.4            metapod_1.8.0             
## [187] dqrng_0.3.0                survMisc_0.5.6            
## [189] R.utils_2.12.2             lazyeval_0.2.2            
## [191] shiny_1.7.4                htmltools_0.5.5           
## [193] KMsurv_0.1-5               sctransform_0.3.5         
## [195] rappdirs_0.3.3             ensembldb_2.24.0          
## [197] glue_1.6.2                 XVector_0.40.0            
## [199] RCurl_1.98-1.12            jpeg_0.1-10               
## [201] gridExtra_2.3              AUCell_1.22.0             
## [203] boot_1.3-28.1              igraph_1.5.0              
## [205] R6_2.5.1                   gplots_3.1.3              
## [207] labeling_0.4.2             km.ci_0.5-6               
## [209] GenomicFeatures_1.52.1     rngtools_1.5.2            
## [211] cluster_2.1.4              Rhdf5lib_1.22.0           
## [213] nloptr_2.0.3               vipor_0.4.5               
## [215] DelayedArray_0.26.6        tidyselect_1.2.0          
## [217] ProtGenerics_1.32.0        raster_3.6-23             
## [219] ggforce_0.4.1              xml2_1.3.4                
## [221] car_3.1-2                  AnnotationDbi_1.62.2      
## [223] future_1.33.0              rsvd_1.0.5                
## [225] munsell_0.5.0              KernSmooth_2.23-21        
## [227] data.table_1.14.8          htmlwidgets_1.6.2         
## [229] RColorBrewer_1.1-3         biomaRt_2.56.1            
## [231] rlang_1.1.1                spatstat.sparse_3.0-2     
## [233] spatstat.explore_3.2-1     lmerTest_3.1-3            
## [235] fansi_1.0.4                beeswarm_0.4.0